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We study dark solitons near potential and nonlinearity steps and combinations thereof, forming rectangular 
barriers. This setting is relevant to the contexts of atomic Bose-Einstein condensates (where such steps can be 
realized by using proper external fields) and nonlinear optics (for beam propagation near interfaces separating 
optical media of different refractive indices). We use perturbation theory to develop an equivalent particle theory, 
describing the matter-wave or optical soliton dynamics as the motion of a particle in an effective potential. This 
Newtonian dynamical problem provides information for the soliton statics and dynamics, including scenarios of 
reflection, transmission, or quasi-trapping at such steps. The case of multiple such steps and its connection to 
barrier potentials is also touched upon. Our analytical predictions are found to be in very good agreement with 
the corresponding numerical results. 
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I. INTRODUCTION 

The interaction of solitons with impurities is a fundamen¬ 
tal problem that has been considered in various branches of 
physics - predominantly in nonlinear wave theory |[I1] and 
solid state physics iH - as well as in applied mathematics 
(see, e.g., recent work ^ and references therein). Especially, 
in the framework of the nonlinear Schrodinger (NLS) equa¬ 
tion, the interaction of bright and dark solitons with ^-like 
impurities has been investigated in many works (see, e.g., 
Refs. iH-ll]). Relevant studies in the physics of atomic Bose- 
Einstein condensates (BECs) have also been performed (see, 
e.g., Refs. IM-B^). as well as in settings involving potential 
wells Q [T^ and barriers in,[13] (see also Ref. (iSD for ear¬ 
lier work in a similar model). In this context, localized impu¬ 
rities can be created as focused far-detuned laser beams, and 
have already been used in experiments involving dark solitons 
|[i 3,[23]- Eurthermore, experimental results on the scattering 
of matter-wave bright solitons on Gaussian barriers in either 
^Lii^ or ^^Rb ll^ BECs have been reported as well. More 
recently, such soliton-defect interactions were also explored 
in the case of multi-component BECs and dark-bright soli¬ 
tons, both in theory ll^ and in an experiment IH]. 

On the other hand, much attention has been paid to BECs 
with spatially modulated interatomic interactions, so-called 
“collisionally inhomogeneous condensates” IH,!!!]; for a re¬ 
view with a particular bend towards periodic such interac¬ 
tions see also Ref. H^j). Relevant studies in this context 
have explored a variety of interesting phenomena: these in¬ 
clude, but are not limited to adiabatic compression of matter- 
waves HI, [13], Bloch oscillations of solitons Hsj], emission 
of atomic solitons llH,!!^], scattering of matter waves through 
barriers iH, emergence of instabilities of solitary waves due 
to periodic variations in the scattering length 132], formation 
of stable condensates exhibiting both attractive and repulsive 
interatomic interactions 1133 ^. solitons in combined linear and 
nonlinear potentials (UllSD, generation of solitons and 


vortex rings ii , control of Earaday waves ll^ . vortex dipole 
dynamics in spinor BECs S, and others. 

Here, we consider a combination of the above settings, 
namely we consider a one-dimensional (ID) setting involv¬ 
ing potential and nonlinearity steps, as well as pertinent rect¬ 
angular barriers, and study statics, dynamics and scattering 
of dark solitons. In the BEC context, recent experiments 
have demonstrated robust dark solitons in the quasi-ID set¬ 
ting iIt^] . In addition, potential steps in BECs can be realized 
by trapping potentials featuring piece-wise constant profiles 
(see, e.g., Refs. HI, SI] and discussion in the next Section). 
Eurthermore, nonlinearity steps can be realized too, upon em¬ 
ploying magnetically |[46] or optically induced Eeshbach 
resonances, that can be used to properly tune the interatomic 
interactions strength - see, e.g., more details in Refs, ll^ l^ 
and discussion in the next Section. 

Such a setting involving potential and nonlinearity steps, 
finds also applications in the context of nonlinear optics. 
There, effectively infinitely long potential and nonlinearity 
steps of constant and finite height, describe interfaces sepa¬ 
rating optical media characterized by different linear and non¬ 
linear refractive indices H^]. In such settings, it has been 
shown iH^-lsl] that the dynamics of self-focused light chan¬ 
nels - in the form of spatial bright solitons - can be effectively 
described by the motion of an equivalent particle in effective 
step-like potentials. This “equivalent particle theory” actu¬ 
ally corresponds to the adiabatic approximation of the per¬ 
turbation theory of solitons while reflection-induced ra¬ 
diation effects can be described at a higher-order approxima¬ 
tion |[53,[5 i 1]. Note that similar studies, but for dark solitons 
in settings involving potential steps and rectangular barriers, 
have also been performed - see, e.^Ref. for an effec¬ 
tive particle theory, and Refs. |[54 - [^ for numerical studies 
of reflection-induced radiation effects. However, to the best 
of our knowledge, the statics and dynamics of dark solitons 
near potential and nonlinearity steps, have not been system¬ 
atically considered so far in the literature, although a special 
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version of such a setting has been touch upon in Ref. ll^ . 

It is our purpose, in this work, to address this problem. In 
particular, our investigation and a description of our presen¬ 
tation is as follows. First, in Sec. II, we provide the descrip¬ 
tion and modeling of the problem; although this is done in 
the context of atomic BECs, our model can straightforwardly 
be used for similar considerations in the context of optics, as 
mentioned above. In the same Section, we apply perturbation 
theory for dark solitons to show that, in the adiabatic approx¬ 
imation, soliton dynamics is described by the motion of an 
equivalent particle in an effective potential. The latter has a 
tanh-profile, but - in the presence of the nonlinearity step - 
can also exhibit an elliptic and a hyperbolic fixed point. We 
show that stationary soliton states do exist at the fixed points 
of the effective potential, but are unstable (albeit in different 
ways, as is explained below) according to a Bogoliubov-de 
Gennes (BdG) analysis |[^|58|] that we perform; we also use 
an analytical approximation to derive the unstable eigenvalues 
as functions of the magnitudes of the potential/nonlinearity 
steps. In Sec. Ill we study the soliton dynamics for various 
parameter values, pertaining to different forms of the effective 
potential, including the case of rectangular barriers formed by 
combination of adjacent potential and nonlinearity steps. Our 
numerical results - in both statics and dynamics - are found 
to be in very good agreement with the analytical predictions. 
We also investigate the possibility of soliton trapping in the 
vicinity of the hyperbolic fixed point of the effective poten¬ 
tial; note that such states could be characterized as “surface 
dark solitons”, as they are formed at linear/nonlinear inter¬ 
faces separating different optical or atomic media. We show 
that quasi-trapping of solitons is possible, in the case where 
nonlinearity steps are present; the pertinent (finite) trapping 
time is found to be of the order of several hundreds of millisec¬ 
onds, which suggests that such soliton quasi-trapping could 
be observable in real BEG experiments. Finally, in Sec. IV we 
summarize our findings, discuss our conclusions, and provide 
provide perspectives for future studies. 


II. MODEL AND ANALYTICAL CONSIDERATIONS 


A. Setup 


As noted in the Introduction, our formulation originates 
from the context of atomic BECs in the mean-field picture 
m. We thus consider a quasi-ID setting whereby matter 
waves, described by the macroscopic wave function T^(x, t), 
are oriented along the x-direction and are confined in a 
strongly anisotropic (quasi-ID) trap. The latter, has the form 
of a rectangular box of lengths '> Ly = Lz = L±, with 
the transverse length L_\_ being on the order of the healing 
length Such a box-like trapping potential, 14 (x), can be 
approximated by a super-Gaussian function, of the form: 


Vb{x) = Fo 


1 — exp 




( 1 ) 


where Vq and w denote the trap amplitude and width, respec¬ 
tively. The particular value of the exponent 7 ^ 1 is not 


especially important; here we use 7 = 50. In this setting, our 
aim is to consider dark solitons near potential and nonlinearity 
steps, located at x = L. To model such a situation, we start 
from the Gross-Pitaevskii (GP) equation 


'"m = 


- + + V{x) 




( 2 ) 


Here, T^(x, f) is the mean-field wave function, m is the atomic 
mass, V(x) represents the external potential, while giD{x) = 
{9/4:L\)gsD is the effectively ID interaction strength, with 
gsD = 4:7rh^a{x)/m being its 3D counterpart and a{x) being 
the scattering length (assumed to be a > 0, Vx, corresponding 
to repulsive interatomic interactions). The external potential 
and the scattering length are then taken to be of the form: 


V{x) 

= 14 W+ 

X < L 
x>L' 

( 3 ) 

a{x) 

_ Jq;l, x<L 
Rr, x>L' 


( 4 ) 


where Vl,r and aL,R are constant values of the potential and 
scattering length, to the left and right of x = L, where respec¬ 
tive steps take place. 

Notice that such potential steps may be realized in present 
BEG experiments upon employing a detuned laser beam 
shined over a razor edge to make a sharp barrier, with the 
diffraction-limited fall-off of the laser intensity being smaller 
than the healing length of the condensate; in such a situa¬ 
tion, the potential can be effectively described by a step func¬ 
tion. On the other hand, the implementation of nonlinear¬ 
ity steps can be based on the interaction tunability of spe¬ 
cific atomic species by applying external magnetic or optical 
fields 1^ . For instance, confining ultracold atoms in an 
elongated trapping potential near the surface of an atom chip 
ll^ allows for appropriate local engineering of the scatter¬ 
ing length to form steps (of varying widths), where the atom- 
surface separation sets a scale for achievable minimum step 
widths. The trapping potential can be formed optically, pos¬ 
sibly also by a suitable combination of optical and magnetic 
fields (see Ref. ll^ for a relevant discussion). 

Measuring the longitudinal coordinate x in units of ^/2^ 
(where i = h/^/2mnglD is the healing length), time t in units 
of a/ 2 (^/Cs (where Cg = \/giD'nlm is the speed of sound and 
n is the peak density), and energy in units of giD'n, we cast 
Eq. © to the following dimensionless form (see Ref. (6^ ): 



1 d^u 

2 dx‘^ 


a{x) 

q^l 


\u\^u -h V{x)u, 


( 5 ) 


where u = ^/n^. Unless stated otherwise, in the simulations 
below we fix the parameter values as follows: Vb = 10 and 
w = 250 (for the box potential), Ul = 0 and Ur = ±0.01 for 
the potential step, as well as aL = 1 and rr G [0.9, 1.1] for 
the nonlinearity step. Nevertheless, our theoretical approach 
is general (and will be kept as such in the exposition that fol¬ 
lows in this section). 

Here we should mention that^. ([5]) can also be applied in 
the context of nonlinear optics 11 ^ : in this case, u represents 
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the complex electric field envelope, t is the propagation dis¬ 
tance and X is the transverse direction, while V (x) and a(x) 
describe the (transverse) spatial profile of the linear and non¬ 
linear parts of the refractive index |[^ . This way, Eq. ([5]) can 
be used for the study of optical beams, carrying dark solitons, 
near interfaces separating different optical media, with (dif¬ 
ferent) defocusing Kerr nonlinearities. 


B. Perturbation theory and equivalent particle picture 

Assuming that, to a first approximation, the box potential 
can be neglected, we consider the dynamics of a dark soliton, 
which is located in the region x < L, and moves to the right, 
towards the potential and nonlinearity steps (similar consider¬ 
ations for a soliton located in the region x > L and moving 
to the left are straightforward). In such a case, we seek for a 
solution of Eq. ([5]) in the form: 


where X = cos0[x — xo{t)] is the soliton coordinate, (p is 
the soliton phase angle {{pl < 7r/2) describing the darkness 
of the soliton, cos p is the soliton depth (0 = 0 and 0 7 ^ 0 
correspond to stationary black solitons and gray solitons, re¬ 
spectively), while xo{t) and dxo/dt = sin 0 denote the soli¬ 
ton center and velocity, respectively. Then, considering an 
adiabatic evolution of the dark soliton, we assume that in the 
presence of the perturbation the dark soliton parameters be¬ 
come slowly-varying unknown functions of time t. Thus, the 
soliton phase angle becomes 0 ^ 0 (f) and, as a result, the 
soliton coordinate becomes X = cos 0(f) (x — xo(f)), with 
dxo{t)/dt = sin 0 (f). 

The evolution of the soliton phase angle can be found by 
means of the evolution of the renormalized soliton energy, 
Eds, given by IHllbl]: 



-ir 


dx. 


( 11 ) 


u{x,t) = Vf/L - exp {-ijxi,t)v{x, f), (6) 

where //l is the chemical potential, and the v{x,t) is the 
wavefunction of the dark soliton. Then, introducing the trans¬ 
formations f ^ (/iL — and X VMl — Vlx, we ex¬ 
press Eq. © as a perturbed NTS equation for the dark soliton: 

Here, the functional perturbation P{v) has the form: 

P{v) = {A^B\v\^)vn{x-L), (8) 


Employing Eq. (fTOl) . it can readily be found that dEds/dt = 
—4 cos^ 0 sin 0 d(j)/dt. On the other hand, using Eq. (|7]) and 
its complex conjugate, yields the evolution of the renormal¬ 
ized soliton energy: dEds/dt = — {Pvt + Pvt)dx, 
where bar denotes complex conjugate. Then, the above ex¬ 
pressions for dEds / dt yield the evolution of 0 , namely 

Inserting the perturbation ® into Eq. (O, and performing 
the integration, we obtain the following result: 


where 1-L is the Heaviside step function, and coefficients A, B 
are given by: 


Vr-Vl d _ 

Ml - Vl <^l 


1 . 


(9) 


These coefficients, which set the magnitudes of the potential 
and nonlinearity steps, are assumed to be small. Such a situ¬ 
ation corresponds, e.g., to the case where ml = 1 , = 0 , 

Vr ^ e, and aR/aj, ^ 1 , where 0 < e <C 1 is a formal 
small parameter (this choice will be used in our simulations 
below). In the present work, we assume that the jump from 
left to right is “sharp”, i.e., we do not explore the additional 
possibility of a finite width interface. If such a finite width 
was present but was the same between the linear and nonlin¬ 
ear interface, essentially the formulation below would still be 
applicable, with the Heaviside function above substituted by 
a suitable smoothened variant (e.g. a tanh functional form). 
A more complicated setting deferred for future studies would 
involve the existence of two separate widths in the linear and 
nonlinear step and the length scale competition that that could 
involve. 

Equation (|7]) can be studied analytically upon employing 
perturbation theory for dark solitons (see, e.g.. Refs, llbll - 
lO ): first we note that, in the absence of the perturbation (jS]), 
Eq. d?]) has a dark soliton solution of the form: 


( 10 ) 


dp 

dt 


— ^[A-\- B) sech^ (L — xo) 

-h -5sech^ (L — xo), 

8 


(13) 


where we have considered the case of nearly stationary (black) 
solitons with cos p ^ 1 (and sin p ^ p). Combining Eq. (O 
with the above mentioned equation for the soliton velocity, 
dxo{t)/dt = sin 0 (f), we can readily derive the following 
equation for motion for the soliton center: 


d‘^xo _ dW 
dP dxo ’ 


(14) 


where the effective potential W (xq) is given by: 


13"(xo) 


- (2A + B) tanh (L — xq) 

8 ^ 

—5tanh^ (L — xo). (15) 


C. Forms of the effective potential 

The form of the effective potential suggests that fixed 
points, where - potentially - dark solitons may be trapped, ex¬ 
ist only in the presence of the nonlinearity step (B ^ 0). I.e., 
it is the competition between the linear and nonlinear step that 


v{x^ t) = cos 0 tanh X i sin 0 , 
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FIG. 1: Sketch showing domains of existence of fixed points of the 
effective potential W{xo) (depicted by gray areas) for A > 0 (blue 
line) and A < 0 (red line). The insets I —III (IV —VI) show the form 
of IV(xo), starting from - and ending to - a small finite value of non¬ 
linearity step B, which is gradually decreased (increased) for A > 0 
(A < 0), cf. black arrows. Small rectangular (yellow) points indi¬ 
cate parameter values corresponding to the forms of W (xo) shown 
in the insets I — VI. 


enable the presence of fixed points and associated more com¬ 
plex dynamics; in the presence of solely a linear step, the dark 
soliton encounters solely a step potential, similarly to what is 
the case for its bright sibling ; see also below. 

In fact, in our setting it is straightforward to find that there 
exist two fixed points, located at: 


xo± 


1 (-AtV-B{2A + B) 

2 ( A + B 


(16) 


for B{2A A B) < 0 and —2A < B < —A, for A > 0 , or 
—A < B < —2A, for A < 0 . In Fig. [T] we plot B(2A + B) 
as a function of 5, for ^4 > 0 (blue line) and < 0 (red 
line). The corresponding domains of existence of the fixed 
points, are also depicted by the gray areas. Insets show typical 
profiles of the effective potential W{xq), for different values 
of B, which we discuss in more detail below. From the figure 
(as well as from Eq. ([16]) itself), the saddle-center nature of 
the bifurcation of the two fixed points, which are generated 
concurrently “out of the blue sky” is immediately evident. 

First, we consider the case of the absence of the nonlinear¬ 
ity step, 5 = 0 , as shown in the insets I and IV of Fig. [T] for 
^4 > 0 and ^4 < 0, respectively. In this case, W{xo) assumes 
a step profile, induced by the potential step. This form is pre¬ 
served in the presence of a finite nonlinearity step, 5/0, 
namely for — ^4 < 5 < 0 and 0 < B < —^4, in the cases 
^4 > 0 and ^4 < 0 , respectively. 

A more interesting situation occurs when the nonlinearity 
step further decreases (increases), and takes values —2A< 
B < —A for A > 0 , or —A < B < —2A for A < 0 . In this 
case, the effective potential features a local minimum (maxi¬ 
mum), i.e., an elliptic (hyperbolic) fixed point, in the region 
X < 0 (x > 0 ) for A > 0 emerge (as per the saddle-center 
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FIG. 2: (Color online) Left panel: density profile of the stationary 
soliton (blue line) at the hyperbolic fixed point xo+ = 0.66, as found 
numerically, using the ansatz 'Us(x) = [1 — V(x)]^/^ tanh(x) in 
Eq. ([it]), for an/aL = 0.985, Vr = 0.01, 14 = 0, /xl = 1; green 
line illustrates the corresponding effective potential W{xo). Right 
panel: corresponding spectral plane (ur, (Xi) of the corresponding 
eigenfrequencies, illustrating an exponential growth due to an imag¬ 
inary eigenfrequency pair. 


bifurcation mentioned above) close to the location of the po¬ 
tential and nonlinearity steps, i.e., near x = 0 ; a similar sit¬ 
uation occurs for A < 0 , but the local minimum becoming a 
local maximum, and vice versa. The locations xq± of the fixed 
points are given by Eq. (HSll; as an example, using parameter 
values Vl = 0, Vr = —0.01, aj, = 1 and Q(r = 1.015, we 
find that xo+ = 0.66 (xq- = — 0 . 66 ) for the elliptic (hyper¬ 
bolic) fixed point. 

As the nonlinearity step becomes deeper, the asymptotes 
(for X —> ±oc) of W(xq) become smaller and eventually van¬ 
ish. For fixed 14 = 0 (and /xl = 1), Eq. (fTSl) shows that this 
happens for B = —{3/2) A; in this case, the potential features 
a “spiky” profile, in the vicinity of x = 0 (see, e.g., upper 
panel of Fig. [ 6 ] below). For B < —{3/2)A, the asymptotes 
of W{xo) become finite again, and take a positive (negative) 
value for X < 0 , and a negative (positive) value for x > 0 , 
in the case A > 0 (A < 0). The spiky profile of W{xo) in 
the vicinity of x = 0 is preserved in this case too, but as B 
decreases it eventually disappears, as shown in the insets III 
and VI of Fig.[T] 


D. Solitons at the fixed points of the effective potential 

The above analysis poses an interesting question regard¬ 
ing the existence of stationary solitons of Eq. © at the fixed 
points of the effective potential. To address this question, we 
use the ansatz xx(x, t) = ex.p{—it)vs{x), for a stationary soli¬ 
ton Vs{x), and obtain from Eq. ([5]) the equation: 

1 (Pvs a{x). ,9 , 

■Us = ---^H - -\vs\ 'l’s + V{x)Vs. (17) 

2 dx^ Q(l 

Notice that we have assumed without loss of generality a unit 
frequency solution; the formulation below can be used at will 
for any other frequency. The above equation is then solved 
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FIG. 3: (Color online) Same as Fig.[2l but for a soliton located at the 
elliptic fixed point xo- = —0.66; this state is found using the initial 
ansatz Vs{x) = [1 — V {x)]^^‘^ tanh(a; + 0.2). The spectral plane in 
the right panel illustrates an oscillatory growth due to the presence 
of a complex quartet of eigenfrequencies. 


numerically, by means of Newton’s method, employing the 
ansatz (for L = 0): 

Vs{x) = [1 — V{x)]^^‘^ tanh(x — xq). (18) 

As shown in the left panel of Fig.[2l assuming an ansatz within 
Eq. (fTSl) in which the soliton is initially placed at xq = 0 , we 
find a steady state exactly at the hyperbolic fixed point xo+ = 
0.66, as found from Eq. (O. On the other hand, the left panel 
of Eig. [3] shows a case where the initial guess is assumed in 
Eq. ([TSb to have a soliton positioned at xq = —0 .2, which 
leads to a stationary soliton located exactly at the elliptic fixed 
point xq- = —0.66 predicted by Eq. (HSl). 

It is now relevant to study the stability of these station¬ 
ary soliton states, performing a Bogoliubov-de Gennes (BdG) 
analysis We thus consider small perturbations of 

Vs{x), and seek solutions of Eq. ([T71) of the form: 

u{x,t) = [vs{x) -h (5 {b{x)e-^^^ + c{x)e^^^)] , (19) 

where (6(x), c{x)) are eigenmodes, cc = -h iuji are (gen¬ 
erally complex) eigenfrequencies, and 6 1. Notice that 

the occurrence of a complex eigenfrequency always leads to 
a dynamic instability; thus, a linearly stable configuration is 
tantamount to = 0 (i.e., all eigenfrequencies are real). 

Substituting Eq. ([T9l) into Eq. ([5]), and linearizing with re¬ 
spect to 6, we derive the following BdG equations: 


^-1 + 2 


a(x) 9 
q^l 


^-1 + 2 


a(x) 9 
q^l 


6 -h 

c-h 


a{x) 

OiL 

a{x) 

C^L 


VgC = ujb, 

Vgb = —ccc, 


( 20 ) 

( 21 ) 


where H = — (1/2)9^ -h is the single particle operator. 
This eigenvalue problem is then solved numerically. Exam¬ 
ples of the stationary dark solitons at the fixed points xq± of 
the effective potential W, as well as their corresponding BdG 
spectra, are shown in Eigs. [2] and [3l It is observed that the 
solitons are dynamically unstable, as seen by the presence of 



FIG. 4: (Color online) Top panel: the imaginary part of the eigen¬ 
frequency, cji, as a function of 1 — B (with B < 0), for a soliton 
located at the hyperbolic fixed point, x = xo+. Middle and bottom 
panels show the dependence of imaginary and real parts, Ui and cjr, 
of the eigenfrequency on 1 — B, for a soliton located at the elliptic 
fixed point, x = xq-, i.e., the case that leads to an eigenfrequency 
quartet. Solid blue curves correspond to the analytical prediction [cf. 
Eqs. (|2^ and ([24l)k blue circles depict numerical results, while yel¬ 
low squares depict eigenfrequency values corresponding to the cases 
shown in Figs. [2| and [3] For the top and middle panels A = 0.01, 
while for the bottom panel A = —(2/3)5; in all cases, /xl = 1. 


eigenfrequencies with nonzero imaginary part in the spectra, 
although the mechanisms of instability are distinctly different 
between the two cases (of Figs. [2| and [3]). 

To better understand these instabilities, and also provide an 
analytical estimate for the relevant eigenfrequencies, we may 
follow the analysis of Ref. |[^ ; see also Ref. for applica¬ 
tion of this theory to the case of a periodic, piecewise-constant 
scattering length setting. According to these works, solitons 
persist in the presence of the perturbation P{v) of Eq. ([5]) (of 
strength A, B e) provided that the Melnikov function con- 
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dition 

/ + ^ QP(y) 

— - —^sech^(x — xo)dx = 0, (22) 

possesses at least one root, say xq. Then, the stability of 
the dark soliton solutions at xo± depends on the sign of the 
derivative of the function in Eq. (l22l) . evaluated at xq: an in¬ 
stability occurs, with one imaginary eigenfrequency pair for 
eM"{xo) < 0, and with exactly one complex eigenfrequency 
quartet for eM"{xo) > 0. The instability is dictated by the 
translational eigenvalue, which bifurcates from the origin as 
soon as the perturbation is present. For eM"(xo) < 0, the 
relevant eigenfrequency pair moves along the imaginary axis, 
leading to an immediate instability associated with exponen¬ 
tial growth of a perturbation along the relevant eigendirection. 
On the other hand, for eM"{xo) > 0, the eigenfrequency 
moves along the real axis; then, upon collision with eigenfre- 
quencies of modes of opposite signature than that of the trans¬ 
lation mode, it gives rise to a complex eigenfrequency quartet, 
signaling the presence of an oscillatory instability. The rele¬ 
vant eigenfrequencies can be determined by a quadratic char¬ 
acteristic equation which takes the form 11 ^ . 

A 2 + ^M"{xo) (l - ^) = (23) 

where eigenvalues A are related to eigenfrequencies cc through 
Since the roots of M"{xo) are the two fixed points 
xo±, we may evaluate M"{xo±) explicitly, and obtain: 

M"{xq±) = — 2sech^(xo±) tanh(xo±) 

X [AB tdiuh^{xo±)] . (24) 

To this end, combining Eqs. (|2^ and (l24l) yields an analytical 
prediction for the magnitudes of the relevant eigenfrequen¬ 
cies, for the cases of solitons located at the hyperbolic or the 
elliptic fixed points of W(xq). 

Figure lU shows pertinent analytical results [depicted by 
(red) solid lines], which are compared with corresponding nu¬ 
merical results [depicted by (blue) points]. In particular, the 
top panel of the figure illustrates the dependence of the imag¬ 
inary part of the eigenfrequency (real part of the eigenvalue) 
uji on the parameter 1 — B (with B < 0), for a soliton lo¬ 
cated at the hyperbolic fixed point, x = xo+; this case is 
associated with the scenario M"{xo) < 0 . The middle and 
bottom panels of the figure shows the dependence of uJi and 
ujr on 1 — B, but for a soliton located at the elliptic fixed 
point, X = Xq-; in this case, M"{xq) > 0, corresponding to 
an oscillatory instability as mentioned above. It is readily ob¬ 
served that the agreement between the theoretical prediction 
of Eqs. (1^ and (l24l) and the numerical result is very good; 
especially, for values oil — B close to unity, i.e., in the case 
\B\ < 0.15 where perturbation theory is more accurate, the 
agreement is excellent. 

We should also remark here that a similarly good agreement 
between analytical and numerical results was also found (re¬ 
sults not shown here) upon using as an independent parameter 
the strength of the potential step (^ ^4), instead of the strength 
of the nonlinearity step (^ B), as in the case of Fig. [H 


III. DARK SOLITONS DYNAMICS 

We now turn our attention to the dynamics of dark solitons 
near the potential and nonlinearity steps. We will use, as a 
guideline, the analytical results presented in the previous sec¬ 
tion, and particularly the form of the effective potential. Our 
aim is to study the scattering of a dark soliton, initially lo¬ 
cated in the region x < L and moving to the right, at the 
potential and nonlinearity steps (similar considerations, for a 
soliton located in the region x > L and moving to the left, are 
straightforward, hence only limited examples of the latter type 
will be presented). We will consider the scattering process in 
the presence of: (a) a single potential step, (b) a potential and 
nonlinearity step, and (c) two potential and nonlinearity steps. 

Attention will be paid to possible trapping of the soliton 
in the vicinity of the location of the potential and nonlinear¬ 
ity steps, and particularly at the hyperbolic fixed point (when 
present) of the effective potential. Notice that in the context of 
optics such a soliton trapping effect could be viewed as a for¬ 
mation of surface dark solitons at the interfaces between opti¬ 
cal media exhibiting different linear refractive indices and dif¬ 
ferent defocusing Kerr nonlinearities (or atomic media bear¬ 
ing different linear potential and interparticle interaction prop¬ 
erties at the two sides of the interface). 

A. A single potential step 

Our first scattering “experiment” refers to the case of a po¬ 
tential step only, corresponding to A > 0 and B = 0 (cf. in¬ 
set I in Fig.[T]). In this case, the effective potential has typically 
the form shown in the top panel of Fig. [51 while the associated 
phase-plane is shown in the middle panel of the same figure. 
Clearly, according to the particle picture for the soliton of the 
previous section, a dark soliton incident from the left towards 
the potential step can either be refiected or transmitted: if the 
soliton has a velocity v = dxo / dt, and thus a kinetic energy 

K = ^ sin^ (j) « (25) 

smaller (greater) than the effective potential step A IF = 
W (+oo) — W (— 00 ), as shown in the top panel of Fig. [51 then 
it will be refiected (transmitted). Notice the approximation 
(sin 0 ^ 0 ) here which is applicable for low speeds/kinetic 
energies. This consideration leads to 0 < 0c or 0 > 0c for 
refiection or transmission, where the critical value 0c of the 
soliton phase angle is given by: 

= v'2AW. (26) 

In the numerical simulations, we found that the thresh¬ 
old between the two cases is quite sharp and is accurately 
predicted by Eq. ([26l) . Indeed, consider the scenario shown 
in Fig. [51 corresponding to parameter values Vy = 0, 
Vr = 0.01, OfR = ay and fiy = 1. In this case, we 
find that AIF = 4.99 x 10 “^, which leads to the critical 
value (for refiection/transmission) of the soliton phase angle 
0c = 9.99 X 10“^. Then, for a soliton initially placed at 
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FIG. 5: (Color online) The case of a single potential step, A = 0.01 
and B = 0, corresponding to 14 = 0, 14 = 0.01, qr = ql, and 
/XL = 1. Top panel (a): effective potential W (xo); shown also is the 
potential difference AW = VF(+(X)) — 14(—oo) = 4.99 x 10“^. 
Middle panel (b): corresponding phase plane; inset shows the initial 
conditions (red squares A and B) for the trajectories corresponding 
to reflection or transmission, while stars and crosses depict respec¬ 
tive PDF results. Bottom panel: contour plots showing the evolution 
of the dark soliton density for the initial conditions depicted in the 
middle panel, i.e., xo = —5 and cj) = 9.6 x 10“^ (left), or 0 = 0.1 
(right); note that, here, 0c = 0.099. Thick (blue) solid lines show 
PDF results, while dashed (white) lines depict ODF results. 


Xo = —5, and for initial velocities corresponding to phase 
angles 0 = 9.6 x 10“^ or0 = 0.1, we observe reflection 
or transmission, respectively. The corresponding soliton tra¬ 
jectories are depicted both in the phase plane (xq, dx^jdt) in 
the middle panel of Fig. [5] and in the space-time contour plots 
showing the evolution of the soliton density in the bottom pan¬ 
els of the same flgure (see trajectories A and B for reflection 
and transmission, respectively). Note that stars and crosses in 
the middle panel correspond to results obtained by direct nu¬ 
merical integration of the partial differential equation (PDF), 
Eq. ([5]), while the (white) dashed lines in the bottom panels 
depict results obtained by the ordinary differential equation 
(ODE), Eq. da. Obviously, the agreement between theoreti¬ 
cal predictions and numerical results is very good. 

Here we should recall that in the case where the nonlinear¬ 
ity step is also present (B ^ 0), and when B > —A (for 
A > 0) or B < —A (for A < 0), the form of the effective po¬ 
tential is similar to the one shown in the top panel of Eig.[5l In 
such cases, corresponding results (not shown here) are qual¬ 
itatively similar to the ones presented above (for 4/0 and 
B = 0); in addition, we have again captured accurately the 
velocity threshold for reflection/transmission. 



FIG. 6: (Color online) Similar to Fig. [S] but for a potential and a 
nonlinearity step, A = 0.01 and B = —0.015, corresponding to 
14 = 0, 14 = 0.01, or/ol = 0.985, and /xl = 1. Top and 
bottom panels show the effective potential W (xo) and the associated 
phase plane, respectively; the potential now features an elliptic and 
a hyperbolic flxed point at xq ^ ±0.65 (cf. vertical dashed lines). 
In the phase plane, initial conditions -marked with red squares- at 
points A (xo = —5, 0 = 0.034), B (xo = —5, 0 = 0.022), C 
(xo = —5, 0 = 0.021) and D (xo = —1.3, 0 = 0.002) lead 
to soliton transmission, quasi-trapping, reflection, and oscillations 
around the elliptic flxed point, respectively; asterisks, crosses and 
stars depict PDF results. The four bottom respective contour plots 
show the evolution of the soliton density; again, thick blue lines and 
white dashed lines depict PDF and ODF results, respectively. 


B. A potential and a nonlinearity step 

Next, we study the case where both a potential and a non¬ 
linearity step are present (i.e., 4,5 / 0), and there exist 
flxed points of the effective potential. One such case that 
we consider in more detail below is the one corresponding to 
4 = 0.01 and B = —0.015 (respective parameter values are 
44 = 0, 44 = 0.01, q^r/q^l = 0.985, and /xl = 1). Note that 
for this choice the effective potential asymptotically vanishes, 
as shown in the top panel of Fig.[6l nevertheless, results qual¬ 
itatively similar to the ones that we present below can also be 
obtained for nonvanishing asymptotics of W (xq). 

The effective potential now features an elliptic and a hyper- 
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bolic fixed point, located at xq ^ ^0.65 respectively. In this 
case too, one can identify an energy threshold AW, now de- 
finedas AVF = W{xo-^)—W{—oo) = VF(xo+), needed to be 
overcome by the soliton kinetic energy in order for the soliton 
to be transmitted (otherwise, i.e., for K < AW, the soliton 
is refiected). Using the above parameter values, we find that 
AW = 2.4 X 10“^ and, hence, according to Eq. ([26l) . the crit¬ 
ical phase angle for transmission/reflection is 0c ~ 0.022. In 
the simulations, we considered a soliton with initial position 
and phase angle xq = —b and 0 = 0.034 > 0c, respectively 
(cf. point A in the phase plane shown in the second panel 
of Fig. [6l), and found that, indeed, the soliton is transmitted 
through the effective potential barrier of strength AW. The 
respective trajectory (starting from point A) is shown in the 
second panel of Fig. [6l Stars along this trajectory, as well as 
contour plot A in the same figure, show PDE results obtained 
from direct numerical integration of Eq. 0; as in the case of 
Fig. [21 the (white) dashed line corresponds to the ODE result. 

To study the possibility of soliton trapping, we have also 
used an initial condition at the stable branch, incoming to¬ 
wards the hyperbolic fixed point, namely xq = —b and 
0 = 0c ^ 0.022 (point B in the second panel of Fig. [ 6 ]). In 
this case, the soliton reaches at the location of the hyperbolic 
fixed point (cf. incoming branch, marked with pluses) and ap¬ 
pears to be trapped at the saddle; however, this trapping occurs 
only for a finite time (for t ^ 600). At the PDE level, this can 
be understood by the the fact that such a configuration (i.e., a 
stationary dark soliton located at the hyperbolic fixed point) is 
unstable, as per the analysis of Sec. II.D. Then, the soliton es¬ 
capes and moves to the region of x > 0, following the trajec¬ 
tory marked with pluses for x > xo+ (here, the pluses depict 
the PDE results). The corresponding contour plot B, in the 
third panel of Fig. | 6 l shows the evolution of the dark soliton 
density. Note that, in this case, the result obtained by the ODE 
(cf. white dashed line) is only accurate up to the escape time, 
as small perturbations within the infinite-dimensional system 
destroy the delicate balance of the unstable fixed point. 

For the same form of the effective potential, we have also 
used initial conditions that lead to soliton reflection. In partic¬ 
ular, we have again used xq = —b and 0 = 0.021 < 0 c, 
as well as an initial soliton location closer to the potential 
and nonlinearity step, namely xq = —1.3, and 0 = 0 . 002 . 
These initial conditions are respectively indicated by the (red) 
squares C and D in the second panel of Fig. | 6 l Relevant tra¬ 
jectories in the phase plane, as well as respective PDE results 
(cf. stars and X marks), can also be found in the same panel, 
while contour plots C and D in the bottom panel of Fig. [ 6 | 
show the evolution of the soliton densities. It can readily be 
observed that for the slightly subcritical value of the phase 
angle (0 = 0 . 021 ), the soliton is again quasi-trapped at the 
hyperbolic fixed point, but for a significantly smaller time (for 
t ^ 150). On the other hand, when the soliton is initially 
located closer to the steps and has a sufficiently small initial 
velocity, it performs oscillations, following the periodic orbit 
shown in the second panel of Fig.[ 6 j 

In all the above cases, we find a very good agreement be¬ 
tween the analytical predictions and the numerical results. 
Similar agreement was also found for other forms of the ef- 



FIG. 7: (Color online) Similar to Fig.jS] for a potential and a nonlin¬ 
earity step, but now for A = 0.01 and B = —0.017, corresponding 
to Vl = 0, Ur = 0.01, or/ol = 0.983, and /xl = 1. The effective 
potential W{xq) (top panel), exhibits an elliptic and a hyperbolic 
fixed point, at xo± = ±0.44 (vertical dashed lines). In the associ¬ 
ated phase plane (second panel) shown are initial conditions, for a 
soliton moving to the right, at points A (xo = — 5, 0 = 0.005) and 
B (xo = — 1, 0 = 0.001), as well as for a soliton moving to the left, 
at points C (xo = 5, 0 = 0.031 > 0c ~ 0.030) and D (xo = 5, 
0 = 0.029 < 0c); in the relevant trajectories, stars, X marks, pluses 
and asterisks, respectively, denote PDE results. Corresponding con¬ 
tour plots for the soliton density are shown in the bottom panels, with 
the dashed white lines depicting ODE results. 


fective potential, as shown, e.g., in the example of Fig. [7] (see 
also inset III of Fig. [T]). For this form of W{xq), parame¬ 
ters A and B dxo A = 0.01 and B = —0.017 (for Vl = 0, 
Ur = 0.01, aR/o^L = 0.983, and /xl = 1), while there exist 
again an elliptic and a hyperbolic fixed point, at xo± = ±0.44 
respectively. In such a situation, if a soliton moves from the 
left towards the potential and nonlinearity steps, and is placed 
sufficiently far from (close to) them - cf. initial condition 
at point A (point B) - then it will be transmitted (perform 
oscillations around xq-). On the other hand, if a soliton is 
initially placed at some xq > xo+ and moves to the left to¬ 
wards the potential and nonlinearity steps, it faces an effec¬ 
tive barrier AW (cf. top panel of Fig. [7]), now defined as 
AW = W (xo+) — lU(±oo). In this case too, choosing an an 
initial condition corresponding to the stable branch, incoming 
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towards xo+, i.e., for the critical phase angle 0c ~ 0.03, it is 
possible and achieve quasi-trapping of the soliton for a finite 
time, of the order oft^ 600. As such a situation was al¬ 
ready discussed above (cf. panel B of Fig. [3), here we present 
results pertaining to the slightly supecritical and subcritical 
cases, namely 0 = 0.031 > 0c and 0 = 0.029 < 0c; cf. 
(red) squares C and D in the second panel, and correspond¬ 
ing contour plots in the bottom panel of Fig. [71 It is readily 
observed that, in the former case, the soliton is initially trans¬ 
mitted through the interface; however, it then follows a trajec¬ 
tory surrounding the homoclinic orbit (see the orbit marked 
with plus symbols, which depicts the PDF results, in the sec¬ 
ond panel of Fig. [7]), and is eventually reflected. In the case 
0 = 0.029 < 0c, the soliton reaches xo+, stays there for a 
time t ^ 180, and eventually is reflected back following the 
trajectory marked with asterisks (see second panel of Fig. [7]). 
In all cases pertaining to this form of W{xo), the agreement 
between the analytical predictions and the numerical results is 
very good as well. 


C. Rectangular barriers 



Our analytical approximation can straightforwardly be ex¬ 
tended to the case of multiple potential and nonlinearity steps. 
Here, we will present results for such a case, where two steps, 
located aX x = —L and x = L, are combined so as to form 
rectangular barriers, in both the linear potential and the non¬ 
linearity of the system. In particular, we consider the follow¬ 
ing profiles for the potential and the scattering length: 


+ 

II 

\x\ > L 
\x\<L' 

(27) 

_ \aR, |a;| > L 
[ C/L, \x\ < L 


(28) 


In such a situation, the effective potential can be found fol¬ 
lowing the lines of the analysis presented in Sec. II.B: taking 
into regard that the perturbation P{v) in Eq. (|7]) has now the 
form: 

P{v) = {A + B\v\‘^)v['H{x + L)-'H{x-L)],i29) 

it is straightforward to find that the relevant effective potential 
is given by: 

W{xo) = - (2A + B) [tanh(I/ — xq) tanh(I/ -h xq)] 

-i- ^5[tanh^(I/ — xq) -h tanh^(I/ + xq)]. (30) 

Typically, i.e., for sufficiently large arbitrary values of L, 
the effective potential is as shown in the top panel of Fig. [8l 
in this example, we used L = 5, while A = 0.01 and 
B = —0.015. It is readily observed that, in this case, asso¬ 
ciated with such a potential and a nonlinearity barrier, is an 
effective potential of the form of a superposition of the ones 
shown in Fig.[6l which are now located at ±5. The associated 


FIG. 8: (Color online) The case of two potential and nonlinearity 
steps forming respective rectangular barriers, for L = 5, A = 0.01 
and B = —0.015, corresponding to Vl = 0, Er = 0.01, or/ol = 
0.985, /iL = 1. Top panel (a): the effective potential W{xo) [cf. 
Eq. ([30l)k featuring elliptic fixed points at the origin and at ±5.66, 
and a pair of hyperbolic fixed points at ±4.34. Middle panel (b): the 
associated phase plane; (red) squares A and B depict different ini¬ 
tial conditions, corresponding to quasi-trapping or oscillations, while 
stars and crosses depict respective PDE results. Bottom panels: con¬ 
tour plots showing the evolution of the dark soliton density for the 
initial conditions depicted in the middle panel, i.e., xo = —8.6 and 
0 = 2.2 X 10“^ (left), or xo = —3 and 0 = 3 x 10“^ (right); here, 
as before, dashed (white) lines depict ODE results. 


phase plane is shown in the middle panel of Eig.[8l shown also 
are initial conditions corresponding to soliton quasi-trapping, 
or oscillations around the elliptic fixed point at the origin - 
cf. red square points A and B, respectively. The correspond¬ 
ing soliton trajectories are depicted both in the phase plane 
in the middle panel of Fig. [S] and in the space-time contour 
plots showing the evolution of the soliton density in the bot¬ 
tom panels of the same figure. Note that stars and plus sym¬ 
bols in the middle panel correspond to PDE results, obtained 
in the framework of Eq. (O, while the (white) dashed lines in 
the bottom panels depict ODE results, obtained by Eq. ilAh 
for the potential in Eq. (l30l Obviously, once again, agreement 
between theoretical predictions and numerical results is very 
good. 

An interesting situation occurs as L decreases. To better 
illustrate what happens in this case, and also to make con¬ 
nections with earlier work JH], we consider the simpler case 
of 5 = 0 (i.e., the nonlinearity step is absent). Then, as¬ 
suming that A = 5/(21/) (with b being an arbitrary small 
parameter), and in the limit of L ^ 0, the potential step 
takes the form of a delta-like impurity of strength b. In this 
case, the effective potential of Eq. (l3Ql) is reduced to the form 


































10 



X 0 



400 800 


1 


0.5 


0 


t 


FIG. 9: (Color online) The case of two potential and nonlinearity 
steps forming respective rectangular barriers, forL = 0.1, A = 0.1 
and B = —0.13, corresponding to aR/uL = 0.87, Vr = 0.1 and 
/XL = 1. Top panel: the effective potential VF(xo), featuring a hy¬ 
perbolic fixed point at the origin and a pair of elliptic fixed points 
at ±1.38. Middle panel: the associated phase plane; (red) square A 
depicts an initial condition corresponding to quasi-trapping of the 
soliton, while stars depict respective PDF results. Bottom panel: 
contour plot showing the evolution of the dark soliton density for 
the initial condition depicted in the middle panel, i.e., xo = — 5 and 
0 = 5.8 X 10“^; as before, dashed (white) lines depict ODE results. 


W{xo) = (6/4Wcli^(xo). This result recovers the one re¬ 
ported in Ref. Jsl] (see also Refs, [[sl [13]), where the inter¬ 
action of dark solitons with localized impurities was studied; 
cf. Eq. (16) of that work, but in the absence of the trapping 
potential f/tr- 

In the same limiting case of small L, and for B 0, the 
effective potential has typically the form shown in the top 
panel of Fig. O here, we use L = 0.1, while A = 0.1 and 
B = —0.13, corresponding to ur/ul = 0.87, Vr = 0.1 
and /iL = 1- Comparing this form of W{xo) with the one 
shown in Fig. [H it becomes clear that as L 0, the indi¬ 
vidual parts of the effective potential of Fig. [8] pertaining to 
the two potential/nonlinearity steps move towards the origin. 
There, they merge at the location of the “central” elliptic fixed 
point, which becomes unstable through a pitchfork bifurca¬ 
tion. As a result of this process, an unstable (hyperbolic) fixed 
point emerges at the origin, while the “outer” pair of the ellip¬ 
tic fixed points (cf. Fig. [8]) also drift towards the origin - in 
this case, they are located at ±1.38. 

In the middle panel of Fig.O shown also is the phase plane 
associated to the effective potential of the top panel. As in 
the cases studied in the previous sections, we may investi¬ 
gate possible quasi-trapping of the soliton, using an initial 


condition at the stable branch, incoming towards the hyper¬ 
bolic fixed point at the origin. Indeed, choosing xq = — 5 and 
(j) = 0c = 6.8 X 10“^ (notice that here, the corresponding 
effective barrier A VF = 1.7 x 10“^-cf. top panel of Fig. [9]), 
we find the following: the soliton arrives at the origin, stays 
there for a time t ^ 600, and then it is transmitted through 
the region x > 0. In fact, the corresponding trajectory found 
at the PDE level is depicted by stars in the middle panel of 
Fig. O while the relevant contour plot showing the evolution 
of the soliton density is shown in the bottom panel of the same 
figure. Notice, again, the fairly good agreement between nu¬ 
merical and analytical results. 

We note that for the same parameter values, but for 5 = 0, 
elliptic fixed points do not exist, and the effective potential 
has simply the form of a sech^ barrier, as mentioned above 
(see also work of Ref. Jsl]). In this case, starting from the 
same initial position, xq = —5, and for 0 = 0.1 (correspond¬ 
ing to 0c = V2AVF ^ 0.1), we find that the trapping time is 
t ^ 320, i.e., almost half of the one that was when the nonlin¬ 
earity steps are present (results not shown here). This observa¬ 
tion, along with the results presented in the previous sections, 
indicate that nonlinearity steps/barriers are necessary either to 
facilitate or enhance soliton trapping in such inhomogeneous 
settings. 


IV. DISCUSSION AND CONCLUSIONS 

We have studied matter-wave dark solitons near linear po¬ 
tential and nonlinearity steps, superimposed on a box-like po¬ 
tential that was assumed to confine the atomic Bose-Einstein 
condensate. The formulation of the problem finds a direct 
application in the context of nonlinear optics: the pertinent 
model can be used to describe the evolution of beams, car¬ 
rying dark solitons, near interfaces separating optical media 
with different linear refractive indices and different defocus- 
ing Kerr nonlinearities. 

Assuming that the potential/nonlinearity steps were small, 
we employed perturbation theory for dark solitons to show 
that, in the adiabatic approximation, solitons behave as equiv¬ 
alent particles moving in the presence of an effective poten¬ 
tial. The latter was found to exhibit various forms, ranging 
from simple tanh-shaped steps - for a spatially homogeneous 
scattering length (or same Kerr nonlinearity, in the context of 
optics) - to more complex forms, featuring hyperbolic and el¬ 
liptic fixed points - in the presence of steps in the scattering 
length (different Kerr nonlinearities). 

In the latter case, we found that stationary soliton states do 
exist at the fixed points of the effective potential. Using a 
Bogoliubov-de Gennes (BdG) analysis, we showed that these 
states are unstable: dark solitons at the hyperbolic fixed points 
have a pair of unstable real eigenvalues, while those at the el¬ 
liptic fixed points have a complex eigenfrequency quartet, dic¬ 
tating a purely exponential or an oscillatory instability, respec¬ 
tively. We also used an analytical approximation to determine 
the real and imaginary parts of the relevant eigenfrequencies 
as functions of the nonlinearity step strength. The analytical 
predictions were found to be in good agreement with corre- 
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spending numerical findings obtained in the framework of the 
BdG analysis. 

We then studied systematically soliton dynamics, for a vari¬ 
ety of parameter values corresponding to all possible forms of 
the effective potential. Adopting the aforementioned equiv¬ 
alent particle picture, we found analytically necessary con¬ 
ditions for soliton reflection at, or transmission through the 
potential and nonlinearity steps: these correspond to initial 
soliton velocities smaller or greater to the energy of the effec¬ 
tive steps/barriers predicted by the perturbation theory and the 
equivalent particle picture. 

We also investigated the possibility of soliton (quasi-) trap¬ 
ping, for initial conditions corresponding to the incoming, sta¬ 
ble manifolds of the hyperbolic fixed points (which exist only 
for inhomogeneous nonlinearities). In the context of optics, 
such a trapping can be regarded as the formation of surface 
dark solitons at the interface between dielectrics of different 
refractive indices. We found that trapping is possible, but only 
for a finite time. This effect can be understood by the fact that 
stationary solitons at the hyperbolic fixed points are unstable, 
as was corroborated by the BdG analysis. Thus small pertur¬ 
bations (at the PDE level) eventually cause the departure of 
the solitary wave from the relevant fixed points. Nevertheless, 
it should be pointed out that the time of soliton quasi-trapping 
was found to be of the order of fiOOv^^/c^ in physical units; 
thus, typically, for a healing length ^ of the order of a mi¬ 
cron and a speed of sound Cs of the order of a millimeter- 
per-second, trapping time may be of the order of ^ 850 ms. 
This indicates that such a soliton quasi-trapping effect may 
be observed in real experiments. Note that in all scenarios 
(refiection, transmission, quasi-trapping) our analytical pre¬ 
dictions were found to be in very good agreement with direct 
numerical simulations in the framework of the original Gross- 
Pitaevskii model. 

We have also extended our considerations to study cases 
involving two potential and nonlinearity steps, that are com¬ 
bined so as to form corresponding rectangular barriers. Re¬ 
fiection, transmission and quasi-trapping of solitons in such 
cases were studied too, again with a very good agreement be¬ 
tween analytical and numerical results. In this setting, special 
attention was paid to the limiting case of infinitesimally small 
distance between the adjacent potential/nonlinearity steps that 
form the barriers. In this case, we found that, due to a pitch- 
fork bifurcation, the stability of the fixed point of the effective 
potential at the barrier center changes: out of two hyperbolic 
and one elliptic fixed point, a hyperbolic fixed point emerges, 
and the potential rectangular barrier is reduced to a delta-like 


impurity. The latter is described by a sech^ effective potential, 
in accordance with the analysis of earlier works lISl-B^ . 

Our methodology and results pose a number of interesting 
questions for future studies. First, it would be interesting to 
investigate how our perturbative results change as the poten¬ 
tial/nonlinearity steps or barriers become larger, and/or attain 
more realistic shapes (including steps bearing finite widths, as 
well as Gaussian barriers - cf., e.g., recent work of Ref. IIJ]). 
In the same context, a systematic numerical - and, possibly, 
also analytical - study of the radiation of solitons durin g re - 
fiection or transmission (along the lines, e.g., of Ref. [SJ]) 
should also provide a more complete picture in this prob¬ 
lem. Furthermore, a systematic study of settings involving 
multiple such steps/barriers, and an investigation of the pos¬ 
sibility of soliton trapping therein, would be particularly rele¬ 
vant. In such settings, investigation of the dynamics of mov¬ 
ing steps/barriers could find direct applications to fundamen¬ 
tal studies relevant, e.g., to superfiuidity (see, for instance. 
Ref. ifT^ l. transport of BECs and even Hawking radia¬ 
tion in analog black hole lasers implemented with BECs |[^ . 
Finally, extension of our analysis to higher-dimensional set¬ 
tings, would also be particularly challenging: first, in order 
to investigate transverse excitation effects that are not cap¬ 
tured within the quasi-ID setting, and second to study simi¬ 
lar problems with vortices and other vortex structures. See, 
e.g.. Ref. l[^ for a summary of relevant studies in higher¬ 
dimensional settings, and Ref. ll^ for a recent example of 
manipulation/control of vortex patterns and their formation 
via Gaussian barriers, motivated by experimentally accessible 
laser beams. 
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